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Abstract 

We compare exact and SU (2)-cluster approximate calculation schemes to determine 
dynamics of the second-harmonic generation model using its reformulation in terms 
of a polynomial Lie àlgebra su p d(2) and related spectral representations of the model 
evolution operator realized in algorithmic forms. It enabled us to implement com- 
puter experiments exhibiting a satisfactory accuracy of the cluster approximations 
in a large range of characteristic model parameters. 
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1 Introduction 

During the last half-century a great attention is being paid to studies of differ- 
ent quantum models of nonlinear òptics since they enable to reveal new phys- 
ical effects and phenomena (see, e.g., [1-13] and references therein). In view 
of the Hamiltonian nonlinearity these models are mainly analyzed with the 
help of the numerical calculations [9,19,20] or some linearization procedures 
[4,7,8] which are not adapted to reveal many peculiarities of model dynamics 
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[2,5,9,11-13]. However, recently a new universal Lie-algebraic approach, es- 
sentially improving both analytical and numerical solutions of physical prob- 
lems, has been suggested in [14] and developed in [6,11,12,15] for the class 
of nonlinear quantum models whose Hamiltonians H have invariance groups 
Gi : [Gi, H] = 0. This approach is based on reformulations of models under 
study in terms of dual algebraic pairs (Gi,g D ) [12,14] where dynamic symme- 
try algebras g D are polynomial Lie algebras (PLA) g p d : g D = g p d completely 
describing model dynamics. Specifically, for a wide class of quantum-optical 
models with g D = g pd = su pd {2) both exact and approximate methods were de- 
veloped in [6,11,12,14,15] to obtain adequate (for physical applications) forms 
of the model evolution operators Unit). Exact methods include appropriate 
Lie-algebraic path-integral schemes [12] and s-u p( i(2)-techniques [6,14,15] based 
on solving difference and difference-differential equations which are fairly com- 
plicated for practical calculations because they yield only algorithms rather 
than closed analytical expressions for solving physical problems. Therefore, 
for this aim two analytical approximations of Un{t) were found in [11]. They 
enable to examine model dynamics at cluster quasiclassical levels described 
by "large" vàlues of certain characteristic parameters [11,12]. However, up to 
now systematic examinations of their accuracy and efficiency were absent. 

In the present work we cancel in part this lacuna for the simplest and most 
widely examined (see, e.g. [1,3,4,6-8,16-20] and references therein) example of 
such models describing process of second-harmonic generation which is given 
by the Hamiltonian 



H = h 



üJ a^a + u-íü^a-í + gaf 2 a + ga\a% , (1) 



and by the Fock space L F (2) = Span{ | {ni, n }) oc (af) ni (aQ ) n °|0)} as the 
model Hilbert space L(H); here ciij circ operators of field modes with fre- 
qüències Ui, g are coupling constants. Classical solutions of the model (1) were 
first found (in terms of elliptical functions) in [1], and later its different quan- 
tum features were examined by many investigators (see, e.g., [3,4,6-8,11-20] 
and references therein). Specifically, in [11,17,19] some peculiarities, like the 
collapse-revival phenomena in the Janés- Cummings model [2], were revealed 
for its model dynamics at the quantum and quasiclassical levels. 

The work is organized as follows. After some preliminaries (Section 2) we give 
exact and approximate calculation schemes and routines to determine "diag- 
onal" (spectral) representations of the model evolution operator by means of 
the sií p d(2)-techniques [14,15] (Section 3). In Section 4 we implement com- 
puter experiments to exhibit their efficiency and to investigate the accuracy 
of appropriate quasiclassical approximations obtained in [11]. In Conclusion 
we briefly discuss some of ways to improve the results obtained. 
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2 The su p d(2) formulation of the model and a related "diagonal" 
representation of its dynamics 

As was shown in [12,14], the invariance of the model Hamiltonian (1) with 
respect to the group Gi = C2 x exp(iA-Ri), C2 = {exp(iriVi) : r = 0, 1} 
enables to re-formulate the model under study in terms of two sets of collective 
operators: integrals of motion i?i = (JVi + 2N )/3, K = JVi - 2[iVi/2] (Ni = 
afüi, [A] is the entire part of A) and Gj-invariant dynamic variables V a= o t ±: 

V + = af 2 a , V_ = a+al V = (N, - N )/3, (2) 
They obey the commutation relations 

[V ,V±] = ±V±, [V a , Rx] = = [V a ,K], 

[V., V + ] = <j)(V ; = V(Vb + 1; Ri) - tl>(y ; Ri), 

4>(V ; Ri) = (Ri + 2V )(R 1 + 2V -1)(R 1 + 1- V ), (3) 

that identifies V a as generators of PLA su pd (2) having the Casimir operator 

*(i?o = {K- RJ/2) = ^j(V ; R 1 ) - V + V_ = vp(V + 1; R ± ) - V_V + , 
[V a ^(R )]=0, (4) 

acting on L(H) complementarily to Gi (in view of the relationship K = 2R + 
Ri) and, hence, forming dynamic symmetry àlgebra g D in the dual algebraic 
pair (Gi,g D ) [12,14]. 

In terms of these collective operators the Hamiltonian (1) is expressed in the 
form [14] 

H = h[ày Q + gV+ + g*V- + C(R 1 )\, 

C=(lu 1 + lü )R 1 , A = 2w 1 -u , K,C} = 0, (5) 
and the Hilbert space L(H) is decomposed into the infinite direct sum 

00 

l{h)= £ £^(M, 

fe=0,l s=0 

L(k, s) = Span{|A;, s; f) : i? í=0 ,i s; /) = U \k, s; /)}, (6) 

of ® SM pc i(2)-irreducible s + 1 - dimensional subspaces L(k, s) specified by 
eigenvalues li of the Gi <g> su p d{2) invariant operators Ri. Herewith li are ex- 
pressed through the numbers s, k : lo — (k — s)/3, /1 = (/c + 2s)/3 determin- 
ing, respectively, a maximal population of the fundamental (pump) mode and 
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a minimal population of the harmònic within a fixed "optical atom" L(k,s) 
as it follows from the structure of the new (collective) basis in L(H) [11]: 



\k, s; /) = N(f; k, s)V+\k, s) = |{m = k + 2f, n = s - /}>, 
N(f;k,s) = [s\M/(k + 2f)\(s-f)\] 1 / 2 , 
V \k,s;f) = (lo + f)\k,s;f), 

\k, s) — \ni — k, no — s), V-\k,s) — 0, k — 0,1, s > 0. 



(7) 



Evidently, Eqs. (7) explicitly manifest the su pd (2) cluster structure of the Fock 
states \ni,n ) and specify R as the "lowest" weight operator and \k, s) as the 
"lowest" weight state. 

This su pd (2) cluster reformulation of the model enables one to use the su pd (2) 
formalism for getting representations of the model evolution operator Unit) = 
exp(—iHt/h) which facilitate analysis of the model dynamics and calculations 
of temporal dependences (0(t)) = Tï\UH{t)pU H {t)0\ for arbitrary physical 
observables O; p = p(0) is the density operator of the initial quantum state 
[6,11]. One of such adequate representations has a diagonal form [11,12] 



where E v (k, s) and {\E v (k, s))} are eigenvalues and complete sets of orthonor- 
malized eigenvectors of the Hamiltonian (5) respectively: 



H\E v (k, s)) = E v (k, s) \E v (k, s)), (E v (h, s)\E v ,(k' ', s')) = 5 k ^5 s ^5 vv ,. (9) 



In the basis {\E v (k, s))} quantum expectations (0(t)) are given as follows [11]: 



U H (t) = E e- ü < k ^\E v (k,s)) (E v (k,s)\, 



E v (k, s) = he(k, s; v) = K[C(li) + X v (k, s)], 



(8) 



(0(t)) 



E (E v (k,s)\p\E v ,(k',s'))(E v ,(k', S ')\0\E v (k,s)) x 




(10) 



and, hence, the main problem is in solving the eigenproblem (9). 



4 



3 Lie-algebraic schemes for finding diagonal su p d(2) forms of the 
model evolution operator: an exact algorithm and cluster quasi- 
classical approximations 



In the case, when the PLA structure polynomial ip{Vo; Ri) is quadratic in Vq 
and PLA su p d{2) is reduced to the familiar Lie àlgebra su(2) = Span{Y , Y± : 
[Y ,Y±] = ±1±, [Y_,Y + ] = 2Y" } the eigenvalue problem (9) is solved ex- 
actly with the help of the SU(2) displacement operators Sy(£ = rg/\g\) = 
exp(£F + — £*Y1) [21] in terms of simple analytical expressions [12]. However, 
it is not the case for the model under study in view of the absence of explicit 
expressions for matrix elements (k, s; f\ exp(X)i a>iVi)\k, s; /') [6,15]. Neverthe- 
less, the su pd {2) formalism enables one get convenient (for physical applica- 
tions) calculation schemes, algorithms and analytical expressions for exact and 
approximate solutions of this problem. 

A Lie-algebraic scheme for finding exact solutions of the eigenproblem (9) is 
based on looking for eigenfunctions \E v (k, s)) on each subspace L(k, s) in the 
form [14,15] 



\E v (k, s)) = J2 Q}(k, s)\k, s;f) = J2 Q}(k, s) V(\k, s), 

f=o f=o 

Q}(k,s) = N(f;k,s)Q}(k,s), (11) 

where amplitudes Q v j(k,s) satisfy the orthonormalization and completeness 
conditions: 

Í2Q}(k,s) Q v ;(k, s ) = 5 VV ,, Í2Q}(^ S ) Q v f>( k > s ) = &fr- ( 12 ) 

/=0 v=0 



Then, inserting Eq. (11) for \E v (k, s)) and Eq. (5) for H in Eq. (9) and using 
Eqs. (3), (7), one gets a set of recurrence relations at fixed k — 0, 1; s — 0, 1, . . . 



g*^l + f + l;h)Q} +1 (k, s) = 

= [X v (k, s) - A(f + l )]Q}(k, s)) - gQ^ik, s), 
f,v = 0,...,s, 

i;(l + f + 1; h) = (k + 2/ + 2)(k + 1 + 2/)(s - /). (13) 

These relations along with the boundary conditions Q^_ ± (k, s) = = Q v s+1 {k, s) 
determine amplitudes Q V f(k, s) and eigenenergies E v (k, s) from solutions of the 
Sturm-Liouville spectral problem [6,14] 
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P /+ i(A) = [A - A(/ + k)]Pf(X) - \g\ 2 4>(l + /; Ji)P/-i(A), (14) 
/ = 0,...,s; P (A) = 1, 

P_x(A) = = [A - A(s + Z )]P S (A) - M 2 ^o + s; Zi)P s _i(A), (15) 

for finding non-classical orthogonal (in view of (12)) polynomials 

p (x) = (9*) f QAk,s;\) = (g*)fQ f (k,s;\) 
A ) N*(f;k,s)Q (k,s;\) N(f;k,s)Q (k,s;X) 1 ) 



of the discrete variable À on the non-uniform lattice {X v (k, s)} s v=0 [14]. Indeed, 
Eqs. (14), (15), (16) provide the following algorithm for solving the eigenprob- 
lem (9). 

i) Using the recursive formula (14) with the boundary vàlues from Eq. (15) 
one calculates the polynomial sequence {P/(À)}j =0 . 

ii) Inserting P s _i(A), P S (X) in the last equality in (15) one gets the algebraic 
equation with respect to À; its solution yields the sequence {A„(A;, s) = A^j^Q 
of admissible vàlues of the spectral parameter À and the appropriate energy 
spectrum {E v = E v (k, s)} s v=0 . 

iii) For each value X v (k,s) using {Pj(à)}j =0 and Eq. (16) one finds the se- 
quence {Q V f(k, s) = Qf(k, s; X v )}f =1 of all amplitudes as functions of the only 
undetermined quantity Qo(X v ) = Qo(k, s; X v ) which, in turn, is found from the 
normalization condition of Eqs. (12). 

This algorithm has been realized in [22] with the help of the REDUCE pro- 
cedure SOLVE for s < 100 and the conventional FORTRAN subroutines of 
EISPARK package (with applying the multiprecision package [23]) for s > 100. 

The routine package developed enables us to implement numerical calculations 
of model dynamics for s < 160 but it is unsuitable for practical calculations 
with larger s in view of multiprecision computer limitations. Therefore in [11] 
an approximate analytical solution of the problem (9) has been suggested. It 
is given by the SU (2)-quasiclassical eigenfunctions 

\E?(k, s; $ = rg/\g\)) = exp(£r + - ?Y-)\k, s; v) = £-^(01*, s; /),(i 7 ) 

/ 

Y = V -l -j, Y + = V + [2(s + 2Y + 2k + l)]- 1 i = (¥_) + , 2j = s, 
Y ±fi G su(2), S} v (0=Qr P (k,s) = (^) f ~ v dl J+L _ J+v (2r), (18) 

and eigenenergies 



Ey C (k, s; = (Er(k, s; Í)\H\K{K s; 0> = n[C(h) + Xf(k, s; r)], 
Xf (k, s; r) = A[j + l - (j - v) cos2r] + 
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+2\g\ £ y/( 8 -/)(/+ l)2(2fc + 1 + 2/) dL 
^ A[j + l — (j — cos2r] — 

-2|#|(j - u) sin2r v /2[s + 2fc + 1 + (-s + 2u) cos2r] = A^ m/ (fc, s; r), (19) 

where d J mn (2r) are the SU (2) ei- functions expressed in terms of the Gauss hy- 
pergeometric function 2 F\ [11]. Approximate vàlues À£ m ^(/c, s; r) in (19) are 
calculated in the cluster mean-field approximation: (k, s; f\F(Y a )\[k, s; f) = 
F((k, s; f\Y a \k, s; /)), and vàlues of the parameter r in (17)-(19) are found 
from energy-stationarity-conditions and/or from minimizing a proximity mea- 
sure between exact Hamiltonian H and its SU (2) - quasiclassical approxima- 
tion 

/f c (0 = £ Er(k,s;0 |£T(M;0) <3T(M;0I- (20) 

í),fe,S 

A Standard measure for such estimates on the subspaces L(/c, s) is defined with 
the help of the unitarily invariant euclidean operator norm [24] as follows [1 1] 

A 2 (h Q , _ Tr (Ks) {H - H^)f _ £, [(À,(M) 2 -(Af(M;r)) 2 ] 
M ' j " Tr (M (ff-C7(Z 1 )) 2 " E,(A„(M) 2 ' 1 j 



This approximation has been used in [11] for calculating approximate expres- 
sions (Y (t)) qc of the temporal dependences (Y (t)) determining, in accordance 
with Eqs. (2), (18), the dynamics of the field-mode populations: 

(N (t)) = S --(Y (t)), (N 1 (t))=s + k + 2(Y (t)), (22) 



for different types of initial states; here s = Tr[p(Ri — Ro)], k = Tr[pK] 
and the quantity (Y (t)) is calculated with the help of Eq. (10) for O = 
Y . Specifically, in the case of the 577(2)-quasiclassical cluster initial state 
\)c of the form (17), (18) with £ = 0,v — 0, belonging to a fixed "optical 
atom" L(k,s) with (jV o (0)) = s, (iVi(O)) = k, Eq. (10) yields an approximate 
analytical expression [11] 



(Y (t))% » ~ {1 + (s - l)A kt8 (t) cos[íl L (k, s)t - (s - l)$ M (í)]} (23) 

with 



nL(fc,s) = 4|W(i-i)(| + fc + i), íï,(m) = MglV ° 



sV2s + 4k + 2' 
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A k ,s(t) 



cos 2 fli(k, s)t + 



sin 2 Qi(k, s)t 



tan 



tanfli(k, s)t 



,(24) 



which exhibits a high-frequency s)) periòdic dynamics with a slow 

(Qi(k, s)) periòdic modulation in the phase and amplitude, i.e. an occurrence 
of a specific temporal coherent structure (described in terms of elliptic func- 
tions too) [11,12]. At the same time for general initial states |) G L(H) and 
having non-zero projections on all subspaces L(k,s), e.g., for Glauber coher- 
ent states [11], analogous calculations lead to series containing weighted sums 
of terms like those given by Eqs. (23), (24) that corresponds to occurrences of 
coherence-decoherence phenomena like " collapse-revivals" revealed in [18,20] 
by means of other methods. 



However, according to the general quasiclassicality theory [25] all approxima- 
tions (17)- (19) (and, hence, (23), (24)) are vàlid only for large vàlues of s, and, 
besides, the measure (21) gives only a global rather than local characteristic 
of the approximate energy spectra {E% c (k, s; £)} that does not allow to feel 
their important symmetry properties and local peculiarities related to "energy 
errors" 

AE v (k, s)=h [X v (k, s) - Xf(k, s; r)] = 5E v (k, s) ■ E v (k, s). (25) 



Therefore, we implemented numerical comparisons of both exact and approx- 
imate results in order to estimate the applicability range of the quasiclassical 
approximation (17)-(19). 



4 Comparison of exact and approximate calculations in the reso- 
nance case 



In order to examine the efhciency of calculation schemes and the algorithm 
given above we tested them by means of computer experiments for the res- 
onance case determined by H from Eq. (11) with u = 2uj 1 ,C{Ii) = Ui(k + 
2s), A = 0. 

First of all we calculated exact vàlues Xf(k, s), Q V f(k, s) according to the algo- 
rithm of Section 3 and their approximations Xl m ^(k, s; r,), Sj v (Ç) according to 
Eqs. (17)-(19) for g = 1, k = 0, 1, s = 20, 10 2 , 5 • 10 2 , 10 3 , 10 4 . Vàlues of the fit- 
ting parameter r were determined from energy-stationarity-conditions: 2r\ = 
arccos | (optimizing only the upper part of spectra) [11], 2r 3 = arccosO = | 
(quasi-linear approximation) [14] and from minimizing the proximity measure 
(21): 2r 2 = arccos ("smooth" cluster mean-field approximation) [11]; here- 
with X^(k, s; =Fri) means that we take r = —r\ in the first half of spectra 
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and r = T\ in the second one. To estimate the accuracy of approximations we 
also used non-invariant measures 



Sl( k , s) = ^KK(k, s) - K m f(k, s; r)] 2 



2 



SeJK S) = E» ' (26) 



to characterize more precisely (in comparison with Eq. (21)) energy spectra 
and standard (related to the Fubini-Study mètric [26] in L(k, s)) measures 

cos(S, Q) [k] . v = £ Sj v Q}(k, s), 5 2 ef (k, s;v) = l-\ cos(S, Q) M;í) | 2 (27) 
/ 



(or associated graphic representations via "overlap areas") to estimate an "ap- 
proximation quality" for eigenfunctions. Some of typical results of these nu- 
merical calculations are presented in Table 1 and Figs.1,2. 



Table 1. Multiplets {A„(A; = 0,s)}^ =0 , {A£ m/ (0, s; ri)} s v=0 with the level step 
Av = 10 for s = 10 2 . 



v X v (0,s) 






A^(0,s;r 2 ) 


A^(0,s;r 3 ) 


-1536.9 


-1096.7 


-1545.3 


-1482.4 


-1421.2 


10 -1151.7 


-919.6 


-1205.2 


-1175.2 


-1137.0 


20 -798.1 


-720.0 


-880.0 


-873.3 


-852.7 


30 -480.3 


-499.3 


-570.2 


-576.7 


-568.5 


40 -205.5 


-259.0 


-276.7 


-285.6 


-284.2 


50 0.0 


0.0 


0.0 


0.0 


0.0 


60 205.5 


276.7 


276.7 


280.0 


284.2 


70 480.3 


570.2 


570.2 


554.3 


568.5 


80 798.1 


880.0 


880.0 


822.8 


852.7 


90 1151.7 


1205.2 


1205.2 


1085.5 


1137.0 


100 1536.9 


1545.3 


1545.3 


1342.3 


1421.2 


K = 


10.222 


-12.220 


0.010 


-1.000 


í| = 


2.563 


0.670 


0.806 


0.657 


51 = 


0.670 


0.670 


0.944 


0.657 



As is seen from data given in Table 1 and Fig. 1 we have an acceptable con- 
sent of exact eigenenergies and their approximations (19) at s ^> 1 almost 
everywhere for r = r í=2i3 and r = ^r\. Discrepàncies between exact and ap- 
proximate results in the middle parts of spectra, probably, are due to the 
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Fig. 1. Energy levels X v (k = 0,s)} (sòlid line) and their approximations (19) with 
r = rj = i ; 2,3 plotted against the energy label v for s = 10 2 ( a) ) and s = 10 4 ( b) ). 

availability of the square-root singularities in the model Hamiltonian (5) re- 
written (with the help of Eqs. (18)) in terms of Y a that is, actually, ignored 
in the "smooth" the su(2) - quasiclassical approximation (17)-(19) [12]. (Note 
also that some negative vàlues of 5 2 H are due to using À£ m ^ instead of À-f 
in Eq. (21) and because of calculation errors). However, the approximation 
with r = =pri breaks the orthogonality of eigenfunctions belonging to opposite 
ends of spectra whereas the quasi-linear approximation with r = r 3 leads to 
equidistant spectra within fixed subspaces L(k,s). Therefore, in spite of the 
spectrum symmetry breaking, the most satisfactory quasiclassical approxima- 
tion is given by Eqs. (17)-(19) with r = r 2 that minimizes 5 2 H . Note that the 
spectrum asymmetry at r = r 2 and related shifts between amplitude vàlues 
Q v j(k,s) and Sj v (t;) (see Fig. 2) are due to using "smooth" su(2) - quasiclas- 
sical eigenfunctions (17) [12]. 



Besides the verifications above we also performed calculations of temporal 
dependences of the quantity (Y (t))c and related dynamics of the normal- 
ized average photon numbers (Ni(t))c/s. Herewith exact dependències were 
calculated with the help of routine package above, whereas approximate cal- 
culations were implemented using approximate expressions (17), (18), (19) for 
eigenvalues and eigenfunctions. Results of such calculations for s = 10 2 , k = 
against the dimensionless time r = gt\/2s are plotted in Fig. 3 where we com- 
pare exact results with the quasiclassical approximations obtained with the 
help of Eqs. (17)-(19) with r = r i= i, 2 and (23)-(24). 
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Fig. 2. Amplitudes Q V f{k, s) as functions in / = 0, s at k = 0, v = 0, 50, 100 (sòlid 
lines) and their quasiclassical approximations S^(r2g/\g\) (dashed lines). 
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Fig. 3. Normalized average photon number (No(t))c/s in the fundamental mode 
plotted against time r for k = 0, s = 10 2 : (exact) is results of (22) and their ap- 
proximations (qc) and (r = ^=1,2) calculated by (23) and (17)-(19) with r = Tí = \^- 

Evidently, the graphic representations of Fig. 3 enable us to reveal transpar- 
ently a double-periodic component in the exact multi-frequency dynamics of 
(Y (t)) and (N (t)) that is rather well described by Eqs. (23), (24) or (17)-(19) 
at r = r 2 (and aX r = r\). Note that an availability of this important dynamic 
feature is displayed clearer when the characteristic parameter s increases (in 
accordance with the general quasiclassical theory [25]). 



5 Conclusion 



So, our numerical calculations given in Section 4 show a rather good qualitative 
consent of exact and approximate results at s ^> 1 and at relevant choices of 
the fitting parameter r in (17)-(19). However, partial quantitative discrepàn- 
cies of them require further improvements of the quasiclassical approximations 
used. In particular, the approximate solutions of the eigenproblem (9) can be 
improved by means of: 1) using less smooth (in comparison with (17), (18)) 
generalized coherent states of the su p d(2) àlgebra as quasiclassical eigenfunc- 
tions (cf. [12]) and 2) exploiting the Standard [24] or special (e.g., developed 
in [27,28]) algebraic perturbative and iterative algorithms or modifications of 
the algebraic "dressing" schemes [6]. Then these improvements (along with 
the exact calculation schemes developed above) can be used for a more detail 
analysis (like those implemented in [7-9,16-20]) of the model under consider- 
ation in all ranges of the parameter s and for arbitrary initial states. It is also 
of interest to compare results obtained (and their improvements) with those of 
[11] based on an alternative SU{2) form of Un{{t) as well as with calculations 
performed in [29] using the formalism of the q - deformed Lie àlgebra su q {2). 
The work along these lines is in progress. 
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